% Generate Panel B in Figure 2.
% To generate Panel A in Figure 2, vary "theta" across -0.6:0.2:0.6, compute the local powers, and take the average.

clear
clc

% Designate working directories
dir = '.../Replication package/';
cd(dir);
addpath(genpath(strcat(dir,"Code/src")))
addpath(genpath(strcat(dir,"Datasets")))

% Set seed
seed = 4534546 ;
rng(seed,'twister');

bnd = 5 ;            % bnd > 0 defines the search space as [-bndbnd]^p, where p = dim(gamma).
bR_P = 1000 ;        % bootstrap replication number used in PSO algorithm.
bR_G = 5000 ;        % bootstrap replication number used in grid search method.
sz_Gamma = 200000 ;  % number of grids in grid search method.
splits = 10 ;
lambdaset = [0.5:(-0.1):0.1 0.05 0.0] ;  % grid for optimal lambda search.
theta = 0 ;         % hypothesized value of θ (EIS) for which the local powers are simulated
B = 2 ;             % Specify B/sqrt(n) as local deviations from "theta".
                    % Represent the collection of local alternatives against which local powers are simulated.
swarmsize = 200 ;   % size of swarm
tol = 1e-4 ;
pso_nmax = 100 ;    % maximum number of PSO trials
alpha_level = 0.1;  % confidence level is equal to 1 - alpha_level

% Load the data
data = table2array(readtable('USAY.txt'));
T = size(data, 1) - 3;
Y = data(4:end, 6);
X = data(4:end, 8);
Z = [data(2:(end-2), 3:6)];

datau = Y - mean(Y) - theta*(X-mean(X));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Local power simulated based on PSO algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
J = length(lambdaset);
noncentral_term = -(X-mean(X))*B;
bootstrap_null = zeros(bR_P, J);
bootstrap_alt = zeros(bR_P, J);
zeta = normrnd(0,1,[T,bR_P]);  % multiplier sequence
eps_constant = 1e-2 ;  % |gamma| < eps_constant is numerically regarded as 0.

% Implement the multiplier bootstrap
for i = 1:bR_P
    test_null = pen_Bierens_test(zeta(:,i).*datau,Z,bnd,lambdaset,eps_constant,pso_nmax,'swarmsize',swarmsize,'dm',1);
    bootstrap_null(i, :) = test_null;
    test_alt = pen_Bierens_test_noncentral(zeta(:,i).*datau,Z,bnd,lambdaset,eps_constant,noncentral_term,pso_nmax,'swarmsize',swarmsize,'dm',1);
    bootstrap_alt(i, :) = test_alt;
end
bootstrap_cv = quantile(bootstrap_null, 1-alpha_level);
LocalPower_PSO = mean(bootstrap_alt > bootstrap_cv);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Local power simulated based on grid search
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Caution: Large "sz_Gamma" may result in memory issues.
%          In such cases, increase the "splits" parameter.
% Set seed again
rng(seed,'twister');
out_GS = pen_Bierens_test_GS(datau,Z,bnd,lambdaset,'bR',bR_G,'sz_Gamma',sz_Gamma,'gradient',-X+mean(X),'splits',splits,'B',B,'seed',seed,'dm',1);
LocalPower_GS = squeeze(min(out_GS.local_powers,[],1));
LocalPower_GS = LocalPower_GS(end:(-1):1);

% Generate plot for Panel B in Figure 2
figure('visible', 'off')
yyaxis left
plot(-lambdaset,LocalPower_PSO,'-','LineWidth',1.5)
ylim([0, 1])
xlabel('Minus penalization parameter (-lambda)')
ylabel('Local power using PSO')
title('Simulated Local Power using PSO and GS')
yyaxis right
plot(-lambdaset,LocalPower_GS,'-','LineWidth',1.5)
ylim([0, 1])
ylabel('Local power using GS')
hold on
plot(-lambdaset,alpha_level*ones(length(lambdaset),1),'-.','LineWidth',0.5,'Color',[0.7 0.7 0.7])
dim = [.15 .1 .3 .1];

saveas(gcf,strcat(dir,'Figure_2.jpg'))